The road to catastrophe: stability and collapse in 2D driven particle systems 
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Understanding collective properties of driven particle systems is significant for naturally occurring 
aggregates and because the knowledge gained can be used as building blocks for the design of arti- 
ficial ones. We model self propelling biological or artificial individuals interacting through pairwise 
attractive and repulsive forces. For the first time, we are able to predict stability and morphology 
of organization starting from the shape of the two-body interaction. We present a coherent theory, 
based on fundamental statistical mechanics, for all possible phases of collective motion. 
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The swarming of multi-agent systems Q is a fascinat- 
ing natural phenomenon. The patterns formed by many 
self-assembling iroecies pose a wealth of evolutionary |2( 
and biological _EJ Q questions, as well as structural and 
physical 0,0,0] ones. In more recent years, understand- 
ing the operating principles of natural swarms has also 
turned into a useful tool for the intelligent design and 
control of man-made vehicles 0, . 

One of the main unresolved issues arising both in ar- 
tificially controlled and biological swarms is the ability 
to predict stability with respect to size. If well defined 
spacings amongst individuals exist, swarm size typically 
increases with particle number, in a 'crystal' like fash- 
ion. This is generally true in animal flocks and might 
be a desirable feature in robotic systems. On the other 
hand, natural examples exist of swarms that shrink in size 
as particle number increases. For instance, in the early 
development of the Myxococcus xanthus or Stigmatella 
auriantaca fruiting bodies [lfl ] two-dimensional bacterial 
vortices arise and grow until the vortices collapse inward, 
individual cells occupy the central core and a localized 
three dimensional structure appears. Although many 
swarming systems have been studied, and in some cases 
specific phase transitions have been observed 
systematic prediction of whether a swarm will collapse 
on itself or not as the number of constituents increases, 
has been lacking. 

In this Letter, we apply fundamental principles from 
statistical mechanics to accurately predict the geome- 
try and stability of swarming systems. Specifically, we 
consider N self-propelled particles powered by biological 
or mechanical motors, that experience a frictional force, 
leading to a preferred characteristic speed 01 ■ The par- 
ticles also interact by means of a two-body generalized 
Morse potential. Previous related work 0, showed that 
in some cases localized vortices may form. Here, we ex- 
plore the entire phase space defined by the interaction po- 
tential and predict pattern geometry and stability. The 
identical N particles obey the equations of motion: 
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with the generalized Morse potential given as: 
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Here, 1 < i < N, and £ a , £r represent the range of the at- 
tractive and repulsive part of the potential and C a , C r are 
their respective amplitudes. We remark that our analy- 
sis need not be confined to Morse-type potentials - these 
are simply a mathematical convenience. Combinations 
of other attractive and repulsive potentials lead to col- 
lective behaviors similar to those seen here. In particular, 
a similar analysis can be extended to the interactions of 
Ref. and of Ref . , providing a much more complete 
prediction of stability. 

In the special case of velocity independent forces, i.e. 
a = (3 = 0, Eg ns . IT151 form a typical Hamiltonian system 
with a conserved energy. For a large number of parti- 
cles, according to the ideas of statistical mechanics, the 
behavior of such a system should be described by a finite 
temperature (Maxwell-Boltzmann) distribution with the 
energy in the model determining the temperature [f3 | . A 
Maxwell-Boltzmann description is even more applicable 
when there is some mechanism for exchange of energy 
with the environment. Here, we argue that the param- 
eters a and (3 in the model of Eans. HO provide a local 
mechanism for such energy exchanges. Thus, the gross 
features of our system should be described by a classical 
statistical mechanics model with the temperature param- 
eter determined by a and f3. Of course the detailed local 
behaviors, as well as the large scale collective dynamics, 
may be quite different - and ostensibly more interest- 
ing - than the corresponding system obeying Maxwell- 
Boltzmann statistics. 

In large systems that obey the laws of statistical me- 
chanics it is expected that thermodynamics will emerge 
as size and number of constituents tend to infinity. In 
order to ensure a smooth passage to the thermodynamic 
limit, the microscopic interactions must respect certain 
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FIG. 1: H-Stability phase diagram of the Morse potential. 
The catastrophic regions correspond to parameter ratios £ = 
£ r /ia and C = Cr/Ca for which the thermodynamic limit 
does not exist. Extrema of the potential d m i„ exist only for 
£ > max{l, C} and for I < min{l, C}. In these cases d m ; n = 
l r log (£/C)/(£-l). 

constraints. The most important of these is H-stability: 
for a set of TV interacting particles, the potential energy U 
is said to be H-stable if a constant B > exists such that 
U > —NB This property ensures thermodynamic 

behavior, in particular that the TV particles will not col- 
lapse as TV — > oo. Non H-stable systems are also called 
catastrophic. For Morse type interactions, conditions for 
H-stability are known. For example, if the d-dimcnsional 
integral of the potential is negative, the system is non 
H-stable. In this case, as TV increases, the particles col- 
lapse into a dense body with energy per particle pro- 
portional to TV. On the other hand, for thermodynamic 
systems, the energy per particle will be asymptotically 
constant. The stability phase diagram of the Morse po- 
tential is shown in Fig. 1. As we shall see the H-stability 
or lack thereof is instrumental in predicting the behavior 
of swarming systems obeying the likes of Eons. lTO 

For the full model, we numerically integrate E qns . Hl^l 
using a fourth order Adams-Bashforth method |l4j for 
free boundaries, effectively allowing an infinite range of 
motion for the particles. Initial conditions are chosen 
with localized particles and random velocities. The re- 
sulting behavior is consistent with the stability or catas- 
trophic predictions of Fig. 1 . We discuss system behavior 
in each region of the phase diagram of Fig. 1. 

In the case min{C, £} < 1 the intcrparticle potential is 
of catastrophic nature and globally attractive. For [3 ^ 0, 
particles tend to sustain a constant speed \v%\ 2 ~ a//?, 
while subject to attractive forces. This competition leads 
to non-equilibrium configurational patterns. We dis- 
tinguish three subregions in the min{C, £} < 1 region: 
{£ < C},{£ = C} and {£ > C}, respectively regions I, 
II and III of Fig. 1. In region I a potential minimum 
d m i n exists and the TV particles self-organize by creating 
multi-particle clumps. Within each clump the particles 
travel parallel to each other defining a collective direc- 
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FIG. 2: Catastrophic geometry, (a) Clumps. From left to 
right TV = 40, 100, 150. Clumps coalesce as TV increases. The 
set parameters are a = 1,(3 = 0.5, C' a = £ a = 1, CV = 0.6, £ r = 
0.5. (b) Ring clumping. From left to right N = 40, 100, 200. 
Parameters are the same as in Fig. 2a, but £ r = 1.2. (c) Rings. 
From left to right N = 60, 100, 200. Parameters are the same 
as in Fig. 2a, but C r — 0.5. (d) Ring radius as a function of 
N from numerical data and Eq.0] Parameters are the same 
as in Fig. 2c. Fitting data from Eq.|I]yields R ~ N~ - 52 . 



tion. Because of the rotational velocity, this direction 
changes in time and the clumps rotate about their center 
of mass (Fig. 2a). Catastrophic behavior is evident in the 
fact that as TV increases, the clumped structures shrink 
instead of swelling. Intcrparticle distance also becomes 
smaller and eventually, as TV — > oo clumps lose their co- 
herence and merge. The bisectant {£ = C}, region II, 
is the borderline for the existence of extrema. Here, the 
potential minimum occurs for d m i n = 0, no associated fi- 
nite length scale exists, and rings are developed (Fig. 2c). 
Assuming equidistant particle spacing, the ring radius R 
may be estimated by balancing the centrifugal and cen- 
tripetal forces. An approximate implicit expression for 
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FIG. 3: Snapshots of swarms for different values of N in 
the catastrophic regime defined by region VII of Fig. 1. From 
left to right N = 100,200,300. Note the decrease in the 
vortex area and the dramatic density increase. The chosen 
parameters are: C a = 0.5, CV = l,l a = 2,£ r = 0.5. Self- 
propulsion and friction are fixed at a = 1.6 and (3 — 0.5. 



R is given by: 
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Estimates of R as given by Eq.0] match extremely well 
those obtained numerically as seen in Fig. 2d. Circular 
structures are also seen in the swarms of Ref. 0. For 
{(. > C} in region III of Fig. 1, clumped structures appear 
although there is no minimum in the potential. In par- 
ticular, no intrinsic intcrparticle spacing exists and the 
clumps consist of superimposed particles traveling along 
a ring: this type of collective motion is energetically more 
favorable than uniform spacing among particles. An ex- 
ample of ring clumping is shown in Fig. 2b. 

A clumped ring structure also appears in the {C < 
1 < £} regime of region IV. The observed behavior is 
very similar to what described in the {£ > C} case above, 
with the difference that here the potential defines a max- 
imum, and for low particle numbers the extra constraint 
of avoiding energetically costly intcrparticle spacings has 
to be considered. Region V of Fig. 1 where max{l, C} > 1 
corresponds to the H-stable regime. Here, the interpar- 
ticle potential is characterized by overall repulsive be- 
havior and is minimized by infinite separation. Thus, as 
N — > oo the particles will tend to occupy the entire vol- 
ume at their disposal. The entire region is a 'gaseous' 
phase with particle speed peaked at l^l 2 = a/ (3. 

The most interesting region of the phase diagram is 
defined by {£ < 1 < C} , regions VI and VII of the phase 
diagram. Here, the potential is characterized by short 
range repulsion and long range attraction. A potential 
minimum exists and defines a length scale d m in- The 
CI 2 = 1 curve of Fig. 1 parts the thermodynamically sta- 
ble region VI from the thermodynamically catastrophic 
region VII. Although the main features of the two-body 
potentials are similar, different H-stability properties 
lead to very different self-organizational behaviors in the 
moderate and large particle limits. 

Region VI with \\j\fC < £ < 1} corresponds to ther- 
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FIG. 4: Vortex scalings for the catastrophic Morse po- 
tential. The parameters are set as in Fig. 3. The friction 
term (3 = 0.5. (a) Vortex area as a function of N for 
a = 1.0, 1.6. Note the dramatic decrease with N. (b) Vortex 
area as a function of (a) for various N. From top to bottom 
N = 90, 140, 200, 300, 400, 600. For any fixed a the vortex 
area decreases with N. (c) Inner and (d) Outer radii of the 
catastrophic vortices as a function of a. The particle num- 
bers are the same as in Fig. 4b. Both radii increase with a 
but decrease with N. For large N the inner core disappears. 



modynamic stability. At finite TV, and for various val- 
ues of a/(3, particles approach the characteristic veloc- 
ity \v 2 \ = a /(3 and reach a kinetic energy much greater 
than the confining interaction potential. The swarming 
agents tend to disperse as individuals. For much smaller 
of values of a//3 , the N particles assemble into orga- 
nized structures with well defined spacings, which in the 
large particle limit tend to a finite value. Particles will 
then either swarm coherently in a rigid disk aggregate 
or flock with a finite center of mass velocity, depending 
on the initial conditions. In both cases, the motion is 
rigid body-like and interparticle distances are preserved. 
For a/f3 — * 0, the particles assemble into static, locally 
crystalline structures. 

Region VII where {£ < 1/VC < 1}, corresponds to 
thermodynamic instability; all cases examined in Ref. 
concern this region. As in the previous case, for finite 
N, large values of a/ (3 will lead to a gaseous phase and 
very small values to crystalline structures whose motion 
is rigid body-like. However, quite unlike the H-stable 
scenario discussed above, these structures are unstable 
with respect to particle number, and in the N — > oo 
limit will collapse. At intermediate values of a/ (3, vor- 
tex structures appear with particles traveling close to the 
characteristic speed \vi\ 2 ~ a/(3. Here, vortex size de- 
creases dramatically as a function of particle number as 
seen in Figs. 3,4. Also, for finite N, vortices rotating 
counter-clockwise and clockwise may coexist, depending 
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FIG. 5: Absolute value of the total Morse potential energy, 
|f/toia(|, as a function of TV. Several choices of a are shown 
with fixed /3 = 0.5. The upper curves correspond to the 
catastrophic regime with the same Morse parameters as in 
Fig. 3. The scaling law is —Utotai ~ TV 2,0 for both a — 1 and 
a = 0. The lower curves correspond to the H-stable parame- 
ters C a = 0.5, C r = 1 and £ a = 2,£ r = 1.5. The scaling law is 
— Utotai ~ JV for both a = and a = 0.5. Scaling is linear 
in the H-stable case, and quadratic in the catastrophic one. 



height of the cells and other features not accounted for 
in EansHEfl 

Simple modeling of interacting particles combined with 
thermodynamic reasoning can account for complex and 
subtle phenomena in biological systems. Furthermore, in 
spite of its name, the existence of a catastrophic regime 
might be of great benefit to the development of un- 
manned vehicle technology for the added versatility it of- 
fers. 'Soft core' robot interactions could be implemented 
by specialized cooperations at short distances e.g. robots 
crawling over or rotating about each other or by set- 
ting parameters so that actual individual size is minis- 
cule compared to relevant length scales. For large or 
even moderate number of vehicles the programming of a 
crossover from an H-stable to a catastrophic regime could 
lead the robots to change from a dispersive (searching) 
behavior to convergence at a specified site. 

We thank H. Levine and D. Marthaler for useful dis- 
cussions. We acknowledge support from ARO and NSF 
through grants W911NF-05-1-0112 and DMS-0306167. 



on the initial conditions. In this regime, the occurrence of 
double spiraling is visually most dramatic since it occurs 
within vortices, however double spiraling is a feature of 
the entire catastrophic part of the phase diagram and co- 
existing left and right direction of motions for clumped or 
cquispaccd rings occur as well. Double spirals arc thus a 
strong indication of the non H-stable nature of the poten- 
tial. Another typical feature of the catastrophic regime 
is that energy per particle does not asymptotically reach 
a constant value. This is seen in Fig.[5]where, in the non 
H-stable regime, the total energy scales quadratically, so 
that energy per particle grows linearly. Here, interpar- 
ticle separation (not shown) decreases dramatically as 
N — > oo. For comparison, in the H-stable regime, the to- 
tal potential energy scales linearly with N and energy per 
particle does approach a constant. Likewise, as N — > oo 
interparticle separation is asymptotically constant. 

The observed swarming behavior of M. xanthus is con- 
sistent with the catastrophic regime VII of the Morse po- 
tential where core free vortex structures can arise. The 
absence of a hard component for the interparticle inter- 
action is justified by the fact that M. xanthus cells can 
penetrate each other by crawling. We propose the follow- 
ing qualitative and coarse scenario for the initial stages 
of aggregation. As the number of constituents increases, 
so does particle speed \J a/ (j. This is consistent with the 
observed enhancement of C-signal activity among par- 
ticles that increases motility. The bacterial vortex 
then increases its size with N and double spirals may 
coexist, as reported in the literature 0|- Eventually, 
bacterial speed reaches an upper limit and increasing N 
will lead the vortex to collapsed behav ior, until, finally, 
it evolves into a complex 3D structure |l8( due to finite 
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